package net.crazyadam.tools;

/**
 * modified by <a href="http://crazyadam.net"><em><i>Joseph J.C. Tang</i></em></a> <br/>
 * Email: <a href="mailto:jinntrance@gmail.com">jinntrance@gmail.com</a>
 *
 * @author joseph
 * @version 1.0
 * @since 10/24/12 11:26 AM
 */

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// This file is part of dLife
// Copyright (c) 2010 Grant Braught. All rights reserved.
//
//     dLife is free software: you can redistribute it and/or modify
//     it under the terms of the GNU General Public License as published
//     by the Free Software Foundation, either version 3 of the License,
//     or (at your option) any later version.
//
//     dLife is distributed in the hope that it will be useful,
//     but WITHOUT ANY WARRANTY; without even the implied warranty
//     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//     See the GNU General Public License for more details.
//
//     You should have received a copy of the GNU General Public
//     License along with dLife.
//     If not, see <http://www.gnu.org/licenses/>.
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import dlife.tools.stat.ArraySummaryStatistics;
import java.io.File;
import java.io.FileNotFoundException;
import java.util.ArrayList;
import java.util.Scanner;

/**
 * Computes the normalized cross correlation between two columns of data
 * contained in files. The files may be space or tab delimited. The output will
 * be two tab-delimited columns. The first column indicates the delay of the
 * second series with respect to the first. The second column gives the
 * normalized cross correlation value at the corresponding delay.
 *
 * @author Grant Braught
 * @author Dickinson College
 * @version Sep 10, 2010
 */
public class CrossCorrelation {
    /**
     * Compute the normalized cross correlation between the two columns of data
     * specified and output the resulting series to standard output. This is a
     * slow (time-domain) implementation of the cross correlation. Thus, it is
     * not recommended for large data sets.
     *
     * <p>
     * The autocorrelation is computed as: <code><pre>
     *
     *         Sum( (s1(i)-ms1) * (s2(i-d)-ms2) )
     * cc[d] = -----------------------------------------------
     *         Sqrt( Sum( (s1(i)-ms1)^2 * Sum( s2(i)-ms2)^2) )
     * </pre></code>
     *
     * <p>
     * Where ms1, ms2 are the means of the s1 and s2 series and i ranges from
     * 0...(n-1) with n equal to the number of values in s1 and s2. The value of
     * d ranges from (n-1)...-(n-1).
     *
     * <p>
     * Usage: java dlife.tools.stat.CrossCorrelation f1 c1 f2 c2 delim
     *
     * @param args <p>
     *            args[0]: Name of file containing the first series (s1). <br>
     *            args[1]: Column of data containing the first series. <br>
     *            args[2]: Name of file containing the second series (s2). Can
     *            be the same as args[0].<br>
     *            args[3]: Column of data containing the second series. Can be
     *            the same as args[1] (for autocorrelation).<br>
     * @throws FileNotFoundException if one of the files does not exits.
     * @throws NumberFormatException if args[1] or args[3] cannot be converted
     *             to an integer.
     */
    public static void main(String[] args) throws NumberFormatException, FileNotFoundException {
        if (args.length != 4) {
            printUsage();
        }
        else {
            double[] s1 = readFile(args[0], Integer.parseInt(args[1]));
            double[] s2 = readFile(args[2], Integer.parseInt(args[3]));
            double[] cc = computeCrossCorrelation(s1, s2);

            for (int i = 0; i < cc.length; i++) {
                System.out.println((i - (s1.length - 1)) + "\t" + cc[i]);
            }
        }
    }

    private static void printUsage() {
        System.out
                .println("Usage: java dlife.tools.stat.CrossCorrelation <f1> <c1> <f2> <c2> <delim>");
        System.out.println("     <f1> : Name of file containing the first series. ");
        System.out.println("     <c1> : Column of data containing the first series (0...n)");
        System.out.println("     <f2> : Name of file containing the second series. ");
        System.out.println("     <c2> : Column of data containing the second series (0...n)");
        System.exit(-1);
    }

    private static double[] readFile(String file, int col)
            throws FileNotFoundException {
        ArrayList<Double> list = new ArrayList<Double>();

        Scanner scr = new Scanner(new File(file));
        while (scr.hasNext()) {
            String line = scr.nextLine();
            String[] toks = line.split("\\s"); // white space.
            list.add(Double.parseDouble(toks[col]));
        }

        double[] dvals = new double[list.size()];
        for (int i = 0; i < list.size(); i++) {
            dvals[i] = list.get(i);
        }

        return dvals;
    }

    /**
     * modifed by Joseph,change the method signature private to public
     * @param s1 array 1
     * @param s2 array 2
     * @return the correlation
     */
    public static double[] computeCrossCorrelation(double[] s1, double[] s2) {
        double[] cc = new double[2 * s1.length - 1];

        double meanS1 = ArraySummaryStatistics.getMean(s1);
        double meanS2 = ArraySummaryStatistics.getMean(s2);

        // Compute denominator for normalization.
        double sS1 = 0;
        double sS2 = 0;
        for (int i = 0; i < s1.length; i++) {
            sS1 += Math.pow((s1[i] - meanS1), 2);
            sS2 += Math.pow((s2[i] - meanS2), 2);
        }
        double denominator = Math.sqrt(sS1 * sS2);

        // Compute numerators for cross correlation.
        int n = s1.length;
        int cci = 0;
        // for each delay of s2 w.r.t s1...
        for (int d = (n - 1); d >= -(n - 1); d--) {
            //System.out.println("d: " + d);
            double sS1S2 = 0;
            // for each entry in s1...
            for (int i = 0; i < n; i++) {

                int j = i + d; // corresponding delayed entry in s2.
                //System.out.println("  i: " + i + " j: " + j);

                if (j >= 0 && j < n) {
                    //System.out.println("    "  + (s1[i] - meanS1)*(s2[j] - meanS2));
                    sS1S2 += (s1[i] - meanS1) * (s2[j] - meanS2);
                }
                else {
                    //System.out.println("     0*");
                    sS1S2 += 0;
                }
            }

            cc[cci] = sS1S2 / denominator;
            cci++;
        }

        return cc;
    }
}